%title:             insurance model: parameter space when risk adjustment
%                   is harmful or beneficial
%author:            neale mahoney
%lastest revision:  6/28/14

%% Setup
clear all; clc;

% cd(<folder directory>);

version = 'v3';

%% Primatives
fprintf('Displaying primatives \n\n');

n_i = 16e3;             %number of individuals
n_c = 1e3;              %number of cost shocks per individual
smooth_type = 'lowess'; %moving, lowess, rlowess
bw = .1;                %bandsmooth for moving average smoother

% joint distribution of risk aversion and expected medical costs
%distribution of absoluate risk aversion from Table 3 of Handel, Hendel and Whinston (2013)
m = 4.39e-4;
v = (6.63e-5)^2;
mu_alpha = log((m^2)/sqrt(v+m^2));
sigma_alpha = sqrt(log(v/(m^2)+1));

%gamble interpretation
f = zeros(100,1);
for i = 1:100,
    f(i) = cara_function(0,m,[-100 i],[0 0]);
end

[value index] = min(abs(f));
fprintf('Mean risk adversion makes individual different between 50-50 gamble for (100,%i) and 0 \n', index)

% distribution of health type
m = .9794951;
v = 1.377593^2;
mu_l = log((m^2)/sqrt(v+m^2));
sigma_l = sqrt(log(v/(m^2)+1));

%joint normal
mu = [mu_alpha; mu_l];
rho = 0;
cov = rho*sigma_alpha*sigma_l;
sigma = [sigma_alpha^2 cov; cov sigma_l^2];
X = mvnrnd(mu,sigma,n_i);

%distribution of risk preferences and health type
alpha = exp(X(:,1));
lambda = exp(X(:,2));

%distribution of realized costs
m = 3138.765;
v = 10125.74^2;
mu_c = log((m^2)/sqrt(v+m^2));
sigma_c = sqrt(log(v/(m^2)+1));
rho = .498481;

%distribution of realized cost conditional on health type
c = zeros(n_i,n_c);
part_1 = mu_c;
part_2 = rho*sigma_c*(X(:,2)-mu_l)/sigma_l;

for i = 1:n_i
    part_3 = sqrt((1-rho^2))*sigma_c*randn(1,n_c);
    c(i,:) = exp(part_1 + part_2(i) + part_3);
end

fprintf('Mean absolute risk aversion parameter is %f \n',mean(alpha));
fprintf('Std dev of absolute risk aversion parameter is %f \n',std(alpha));

fprintf('Mean population cost is %i \n',round(mean(reshape(c,n_i*n_c,1))))
fprintf('Std dev of cost is %i \n',round(std(reshape(c,n_i*n_c,1))))

%av 60 plan: low quality plan that covers 60% of medical cost on average
c_oop_60 = min(min(c,500+.4*(c-500)),8e3);  %500 deductible, 8000 out of pocket maximum
c_ins_60 = c - c_oop_60;
c_ins_mean_60 = mean(c_ins_60,2);

fprintf('AV 60 mean plan cost is %i \n',round(mean(c_ins_mean_60)))
fprintf('AV 60 mean oop cost is %i \n',round(mean(reshape(c_oop_60,n_i*n_c,1))))
fprintf('AV 60 AV is %f \n',mean(c_ins_mean_60)/mean(reshape(c,n_i*n_c,1)))

% av 90 plan: high quality plan that covers 90% of medical cost on average
c_oop_90 = min(.1*c,8e3);   %no deductible, 8000 out of pocket maximum
c_ins_90 = c - c_oop_90;
c_ins_mean_90 = mean(c_ins_90,2);

fprintf('AV 90 mean plan cost is %i \n',round(mean(c_ins_mean_90)))
fprintf('AV 90 mean oop cost is %i \n',round(mean(reshape(c_oop_90,n_i*n_c,1))))
fprintf('AV 90 AV is %f \n',mean(c_ins_mean_90)/mean(reshape(c,n_i*n_c,1)))

%allow bankruptcy to cap out of pocket risk
c = min(c,20e3);

%% wtp
fprintf('\n\nEstimating WTP \n\n');

v_grid = 0:10:5e4;
v = zeros(n_i,1);
for i = 1:n_i,

    if mod(i,25) == 0, fprintf('Estimating WTP for observation %i \n',i); end
    v(i) = wtp_function(alpha(i),v_grid,c_oop_90(i,:),c_oop_60(i,:));

end
 
%% Equibriums

theta_grid = 0:.05:1;
wtp_shift_grid = -1e3:50:1e3;

%sigma = 0 is full ra, 1 is no ra, 3 is negative ra;
sigma_grid = [.9 1];

out = zeros(length(theta_grid),length(wtp_shift_grid),length(sigma_grid),7);

for i = 1:length(theta_grid),
    for j = 1:length(wtp_shift_grid),
        for k = 1:length(sigma_grid),
            
            theta = theta_grid(i);
            wtp_shift = wtp_shift_grid(j);
            sigma = sigma_grid(k);
            
            fprintf('Calculating equilibrum for theta %i, wtp shift = %f, sigma = %f \n', theta, wtp_shift, sigma);
            
            %sort by wtp
            X = [v+wtp_shift c_ins_mean_90-c_ins_mean_60];
            X = sortrows(X,-1);
            v_sort = smooth(X(:,1),bw,smooth_type);
            
            mc = smooth(X(:,2),bw,smooth_type);
            q = transpose((1:n_i)./n_i);
            
            out(i,j,k,:) = equilibrium_function_v3(v_sort,mc,q,theta,sigma,bw,smooth_type);
            
        end
        
    end
end


%% Quantity thresholds
q_base          = zeros(length(theta_grid),length(wtp_shift_grid));
ra_excessive    = zeros(size(q_base));
q_ra_excessive  = zeros(size(q_base));

ra_good_firm    = zeros(size(q_base));
q_ra_good_firm  = zeros(size(q_base));

ra_good         = zeros(size(q_base));
q_ra_good       = zeros(size(q_base));

ra_neg          = zeros(size(q_base));
q_ra_neg        = zeros(size(q_base));

for i = 1:length(theta_grid),
    for j = 1:length(wtp_shift_grid),
        
        sigma_base = find(sigma_grid==1);
        q_base(i,j) = out(i,j,sigma_base,1);
        
        cs          = squeeze(out(i,j,:,4));
        ps          = squeeze(out(i,j,:,5));
        transfer    = squeeze(out(i,j,:,6));
        fs          = cs - transfer; 
        ss          = squeeze(out(i,j,:,7));
                
        %some ra good for firm
        %cs - transfer with some ra > cs - transfer at baseline
        ra_good_firm(i,j) = (fs(sigma_base-1) >= fs(sigma_base));
        q_ra_good_firm(i,j) = q_base(i,j)*ra_good_firm(i,j);
        
        %some ra good
        %q some ra > q baseline
        ra_good(i,j) = (out(i,j,sigma_base-1,1) > out(i,j,sigma_base,1));  
        q_ra_good(i,j) = q_base(i,j)*ra_good(i,j);
        
    end
end

ra_good_firm_threshold = max(q_ra_good_firm,[],2);
ra_good_firm_threshold_s = smooth(ra_good_firm_threshold);
ra_good_firm_threshold_s(1) = 1;

ra_good_threshold = max(q_ra_good,[],2);
ra_good_threshold_s = smooth(ra_good_threshold);
ra_good_threshold_s(1) = 1;

%Simulate markets based on the distribution of market power reported in
%Dafny, Duggan and Ramanarayanan (2012)
n_sim = 1000;
x = betarnd(1.665312, 1.256288, n_sim,1);
y = betarnd(19.86286, 46.70169, n_sim,1);

figure(1)
hold all;
scatter(x,y,1,'MarkerEdgeColor',[0.8 0.8 0.8],'linewidth',2);
plot(ra_good_firm_threshold_s,theta_grid,'LineStyle', '-','Color',[0.5, 0.5, 0.5], 'linewidth',2);
plot(ra_good_threshold_s,theta_grid,'LineStyle', '-','Color',[0, 0, 0], 'linewidth',2);
box on

ylabel('Market Power (\theta)','fontsize',12); xlabel('Baseline Quantity (q^*)','fontsize',12);
axis([0 1 0 1])
str = {'Risk Adjustement Reduces', '         Social Surplus'};
text(.6, .8, str ,'FontSize',10)
str = {' Risk Adjustement', 'Reduces Employer +',  ' Employee Surplus'};
text(.2,.7,str,'FontSize',10)
str = {'Risk Adjustement Beneficial'};
text(.1,.1,str,'FontSize',10)

f_out = '../Figures/Risk_adjustment_harmful_v3';
hgsave(f_out)
print('-painters','-dpdf','-r300',f_out);


